clear, close, clc

w = 2*pi;
d = 0.5;

A = [0 1; -w^2 -2*d*w];

dt = 0.1; % time step
T = 10; % time to integrate (duration)

x0 = [2; 0]; % initial condition

% iterate forward Euler scheme
xF(:,1) = x0;
tF(1) = 0;
for i=1:T/dt
    tF(i+1) = i * dt;
    xF(:,i+1) = (eye(2) + A*dt)* xF(:,i);
end
plot(tF,xF(1,:),'k')
xlabel('Time')
ylabel('position')
hold on

% iterate backward Euler scheme
xB(:,1) = x0;
tB(1) = 0;
for k=1:T/dt
    tB(k+1) = k*dt;
    xB(:,k+1) = inv(eye(2)-A*dt) * xB(:,k);
end
plot(tB,xB(1,:),'r')

% try using ode45
[t,y] = ode45(@(t,y) A*y, 0:dt:T,x0);
plot(t,y(:,1),'b')

legend('Forward E','Backward E','ode45')